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Implementations of a model for equilibrium, steady-state ablation boundary conditions 
are tested for the purpose of providing strong coupling with a hypersonic flow solver. The 
objective is to remove correction factors or film cooling approximations that are usually 
applied in coupled implementations of the flow solver and the ablation response. Three test 
cases are considered - the IRV-2, the Galileo probe, and a notional slender, blunted cone 
launched at 10 km/s from the Earth’s surface. A successive substitution is employed and 
the order of succession is varied as a function of surface temperature to obtain converged 
solutions. The implementation is tested on a specified trajectory for the IRV-2 to compute 
shape change under the approximation of steady-state ablation. Issues associated with 
stability of the shape change algorithm caused by explicit time step limits are also discussed. 


I. Introduction 

Exploration missions to planets or moons supporting an atmosphere and return missions to Earth usually 
involve high speed entry at hyperbolic speeds (greater than escape velocity). Atmospheric drag is used to 
slow the vehicle to enable aerocapture or landing and save mass associated with fuel burn of retro-propulsion 
prior to entry. Simulation of the aerothermal environment surrounding a vehicle traveling at hyperbolic 
velocities through a planetary atmosphere requires coupling of the high temperature flow field with the 
thermal protection system (TPS) material response (ablation, shape change) and with the transport of 
energy across the shock layer by radiation. In the case of high mass missions to Mars, atmospheric density is 
relatively thin and supplementary braking concepts, including deployable / inflatable surfaces to increase drag 
or supersonic retro-propulsion need to be engaged, introducing new coupling requirements in the simulation. 

The challenge of coupling is that small perturbations to shared boundary conditions across the gas - solid 
interface can lead to large and destabilizing perturbations in the solution of the multi-physics components. 
Loosely coupled algorithms involving multiple flow solver relaxation steps for every update to the material 
response code were the first approaches to deal with both the complexity and numerical stability in the 
simulations of ablation. The evolution of loosely coupled material response codes to flow solvers has been 
described by Chen and Milos. 1 and Kuntz et. al. 2 The Charring Material Thermal Response and Ablation 
Program 3 (CMA), developed by the Aerotherm Corporation in the 1960s, the Fully Implicit Ablation and 
Thermal Response program 4 (FIAT), developed at NASA Ames Research Center in the 1990’s, and Chaleur, 5 
developed at Sandia and Johnson Space Center in the 2000’s, are one-dimensional material response codes. 
Two-dimensional material response codes, better able to handle recession at sharp corners and nose tips, 
have been developed more recently, including COYOTE 6-8 from Sandia National Laboratories and the Two- 
dimensional Implicit Thermal Response and Ablation (TITAN) Program 1 from NASA Ames. A strongly 
coupled algorithm with material response updates implemented at every update of the flow solver, albeit 
with simplifications in the gas chemistry model, is described by Conti et. al. 9 A more recent strongly coupled 
algorithm, including effects of micromechanical ablation, was presented by Gosse and Candler 10 

Most of these methods require use of a heat transfer coefficient and access to the thermal boundary layer 
edge conditions in order to specify the surface energy balance that connects the flow solver to the material 
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response. 11 See, for example, Eq. (1) in Kuntz et. al. 2 * for coupling to COYOTE or Eqs. (8) and (9) 
in Chen and Milos 1 for coupling to TITAN. These approximations recognize that small perturbations to 
temperature and species concentration at the wall produce large perturbations to heating and diffusion rate. 
Even when line-implicit relaxation is used the abrupt change to the near wall gradients governing heating and 
diffusion can be destabilizing. Use of heat transfer coefficient approximations moderate these perturbations 
and promote convergence of the coupled systems. (Only Gosse and Candler 10 do not indicate the use of 
any approximations in the implementation of a surface energy balance nor do they discuss any algorithm 
challenges to the implementation of the interface condition.) 

As noted above, simulations of aerothermal loads with ablation are evolving to strong coupling, with 
updates to the material response synchronized with updates to the flow solver for every relaxation step. 
Strong coupling provides greater opportunity to relax the solution across the interface and ultimately re- 
move any coefficient approximations. Removing the heat transfer coefficient approximations in the surface 
energy balance enables a more accurate simulation - especially in situations where local geometry includes 
characteristic lengths that are small compared to the boundary layer thickness (e.g. cavities, gap fillers, 
corners). The long term goal in these topologically challenging simulations is to couple multi-dimensional 
material response with locally resolved flow gradients - unencumbered by any constraints to approximate 
heat transfer coefficients. In the process of implementing these strongly coupled algorithm changes stability 
problems were encountered unless very small small relaxation parameters were used that pushed the number 
of relaxation steps required to > 50000 for a typical, axisymmetric simulation. The purpose of this paper 
is to describe how modifications to the boundary condition relaxation enabled more acceptable convergence 
rates. Three cases are presented with large variations in blowing rates, surface temperatures, and convective 
and radiative heating that test the new boundary condition relaxation algorithm. Future tests are planned 
to engage more accurate material response on topologically complex surfaces. 


II. Nomenclature 


Bold face, lowercase variable names refer to vectors. Bold face, uppercase variable names refer to matri- 
ces. Bracketed entry indicates units or quantity used to non-dimensionalize. 


Roman symbols 

Aqi, A02 
Aj 

A R 
J - f n 

C 3 

Ci 

Ci,abl 

Dj 

F 

Fi,j 

Gj,k 

h 

k 

Kj 

m 

Mj 

n 

h 

P 


Jacobian matrices for boundary condition updates 

area of cell face j on surface 

Fourier coefficients 

mass fraction of species j 

elemental mass fraction of element i 

elemental mass fraction of element i in ablation products 
effective diffusion coefficient for species j 

transformation matrix from species continuity to elemental continuity equations 
element of F 

coefficient of In in partial equilibrium relation for non-basis species j 

enthalpy per unit mass 

turbulent kinetic energy 

equilibrium constant for non-base-species j 

blowing rate of ablation products per unit area 

molecular weight of species j 

normal distance to wall 

recession rate 

pressure 
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q 

r 01) r 02 
t 

x,y,z 


heating rate 

vector of residuals in ghost cell 
time 

Cartesian coordinates 


Roman symbols 

T 

V 

Vj 

W0l,W 0 2 


temperature 

velocity normal to surface 
diffusion velocity of species j 
vector of unknowns in ghost cell 


Greek symbols 

a 

e 

£ 

P 

<7 

Xj 


fraction of radiation absorbed at surface 

emissivity 

relaxation factor 

density 

Stefan-Boltzmann constant 
mole fraction of species j 


Subscripts 

0 

1 

c 

C 

cond 

char 

9 

i 

j 

k 

rad 

v 

w 

00 


ghost cell behind boundary 
internal cell adjacent to boundary 
convective 
carbon 

conduction through surface 

property of char 

pyrolized gases 

element index 

species index 

species index 

radiative 

virgin material 

wall 

reference condition in free stream 


Superscripts 

m 

n 


sub-iteration index 
global relaxation index 


III. Conservation Equations 

Program LAURA (flow solver) is used to generate all simulations presented herein. The conservation 
equations are provided in the literature . 12, 13 These include species conservation, momentum conservation, to- 
tal energy conservation, and vibrational-electronic energy conservation if thermal non-equilibrium is modeled, 
and turbulence models. Park’s two-temperature model is employed in the case of thermal non-equilibrium 
to describe modifications to the chemical kinetics. A Free Energy Minimization (FEM) algorithm is used 
in association with elemental species continuity equations in the case of thermochemical equilibrium model- 
ing. In cases where radiative energy transport is included a tangent-slab approximation is used through the 
HARA modules that are fully-coupled within LAURA. 
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IV. Surface Boundary Conditions 


Both pyrolysis and ablation are accommodated in the current analysis. Pyrolysis describes the internal 
decomposition of the solid which releases gaseous species into the shock layer. Ablation describes the 
combination of processes which consume heatshield surface material. 1 An equilibrium surface boundary 
condition requires that all species are equilibrated at the local wall temperature T w with elemental mass 
fractions determined by the elemental continuity equation. Given a system with N s species and N e elements 
the algorithm identifies N e base species, usually atoms, and N s — N e non-base-species which can be created 
from the base set. The transformations from species continuity equations to elemental continuity equations 
and the definition of equilibrium constant Kj(T) and stoichiometric related coefficients Gj t k for each of 
the linearly independent non-base-species j shown below were derived previously. 13 The N s — N e partial 
equilibrium relations for non-base species is written 


species 

Kj(T w ) = ^ Gj t k In pk 

k 


(1) 


where In pk is discretized 0.5(ln pkp + In Pk,i) and where subscript 0 denotes a ghost cell beneath the surface 
and subscript 1 denotes the cell center above the surface in the cell-centered formulation. 

The integrated form of the species continuity equation 13 provides a balance for each element in which the 
convectecl elemental mass flux out of the boundary minus the diffused elemental mass flux into the boundary 
must equal the ablated elemental mass flux specified for the material. It is written 


F it jPjV ~ G,abirn = F itj pD 


dXj 

dn 


Cj J c< 


(2) 


where F t j is a matrix of size (N s ,N e ) that transforms species mass fraction to elemental mass fraction and 
V is velocity normal to the surface. Note that an effective binary diffusion coefficient is applied to model 
diffusive flux. A mass corrected diffusive flux 14 is also applied using J cor = pDj ^ to guarantee that 
the sum of diffusive flux over all species (or over all elements) goes to zero. The term FijpjV is discretized 
0.5(piCj,iVi + FijPjfiVo). Other terms are evaluated as arithmetic averages or as central differences. The 
elemental mass fraction of ablated gas, Ci, a bi, is a property of the thermal protection system. 

In the case of no ablation the sum of Eq. 2 over all elements i is zero. These equations are not linearly 
independent and so the one with the largest elemental mass fraction at the cell center bordering the wall is 
replaced with a specification of normal momentum conservation assuming viscous terms may be ignored. 


A surface energy balance is written 


Qc QQrad T “b Qcond T {Fl char T Fig) fl w FlcharFchar Flghg — 0 


( 4 ) 


In the present formulation assuming steady state ablation with the char front and virgin material front 
receding at the same rate, the elemental mass fractions of the pyrolysis gas, c g y, and the char, c c har,i are 
known so that 


. TTlchar 

'H'char — — 

Pchar 


m 9 

Pv Pchar 


( 5 ) 


from which it follows, using rrii = m c har,i + 


Ci,abl 


( Pchar 
Pv 


Cchar, i ~t" 



Pchar 


C 9,i • 


( 6 ) 


The steady state ablation approximation also assumes that the conduction of energy into the TPS is 
balanced by the energy content in the gas approaching the TPS surface (minus a negligibly small energy 
conduction beyond the pyrolysis front) so that 


Qcond — ITlcfiarhchar T Flghg 


( 7 ) 
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These approximations simplify the surface energy balance Eq. 4 to 

Qc QQrad T 4“ ( J^char 4“ dlg^)H w 0 (3) 

Note that a material response code would be required to solve for q con d and rh g if a steady-state ablation 
approximation is not used. As will be seen subsequently, a steady-state ablation approximation over-predicts 
the recession rate early in the trajectory. 

If there is no ablation then all surface velocities (and blowing rates) are zero. The N s species densities 
are defined with N s — N e equations for non-base-species (Eq. 1), N e — 1 equations for base species through 
elemental conservation (Eq. 2), and one equation for normal momentum conservation (Eq. 3). Surface 
temperature is obtained from the surface energy balance equation (Eq. 8) which, in this case of no ablation, 
is simply an expression of the radiative equilibrium wall temperature approximation. Vibrational - electronic 
temperatures are assumed equal to translational temperature at the surface. These equations define the 
equilibrium surface catalytic boundary condition. 

If there is ablation then an additional boundary condition is needed to define m c har- (Note that m g is 
obtained from the steady-state ablation approximation (Eq. 5, 7) or from a material response code. Also note 
that the char density p c har and the virgin material density p v are assumed known properties of a charring 
ablator.) The closing relation recognizes that the solid surface char participates as a reaction partner with 
the gas. It is implicitly assumed here that a single equilibrium relation between gaseous carbon and solid 
carbon in the char is sufficient to define m c har , regardless of the presence of any other elements in the char. 
More comprehensive models of material response 15 dealing with this (and other) simplifications have not 
been tested here. Therefore, the closing relation to define m c har is expressed 

PC,w = M cKchar,c{T w ) (9) 

A curve fit of the equilibrium constant K c har,c is expressed 

Kahar.c = 1-238 lO 11 ! 1 - 1 - 487 exp ^ ~ 87 ^ 10 - 4 ^ ( 10 ) 

Note that this closing relation has only implicit dependence on m c har ; the blowing rate is adjusted until the 
equilibrium condition is attained at a wall temperature that satisfies the surface energy balance. 

A. Implementation of Boundary Conditions 

It is observed that the mutual sensitivity of dT w , dm, and din pc, o changes with increasing surface temper- 
ature - likely due to diffusion limited presence of oxygen near the wall as blowing rates increase. A surface 
temperature switch T sw is defined to enable convergence. Three conditions with equilibrium surface chem- 
istry are considered: (1) no ablation; (2) ablation with T w < T sw ; and (3) ablation with T w > T sw . T sw is a 
temperature that defines where the implementation algorithm is switched. Numerical experiments indicate 
T sw = 2800K. The boundary conditions at ghost cells are updated in parallel with the interior cells at every 
relaxation step n. As will be shown below, when T w > T sw the expected atomic carbon mass fraction cc at 
the surface starts to rise above a trace level and so the value of cc from the current interior relaxation step 
can be used to update the surface temperature for the next relaxation step. When T w < T sw the expected 
m is small, and the linearized set of boundary conditions, including the implicit dependence of m on the 
species densities, can be solved using a strongly coupled Newton relaxation. 

Note that in this solution procedure the convective heating and diffusion of species to the surface are 
calculated directly from the local gradients. There is no approximation involving film cooling coefficients. The 
drawback of this approach is that small changes in boundary values can lead to large, physically unrealistic 
changes in the convective heating rates and diffusion rates. The off-body points take longer to respond in the 
global solution. The use of line relaxation and relaxation factors on boundary updates enable convergence 
on the order of 20,000 to 40,000 global relaxation steps, including re-gridding and shock alignment for the 
evolving solution. 
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1. No Ablation 


In the case of no ablation the ghost cell unknowns include species densities and temperature. All other 
quantities are explicitly formulated from these unknowns. Temperature is first computed from Eq. 8 as 

T n+1 = £ ' + (1 _ e ^ T u (n) 

with relaxation factor e = 0.001. Densities are computed by the simultaneous solution of Eqs. 1, 2, and 3 
for In pjfi with surface ablation rates and normal velocity set to zero. A sub-iteration over index m of these 
equations with a Newton method is employed to a residual convergence of lrrlO -10 . 

ASlAw 01 = r™ (12) 

Here Aw 0 i = ("« — wjj) is a vector of N s unknowns equal to In pj t0 . r 0 i is the vector of residuals of 
Eqs. 1, 2, and 3 evaluated at T™ +1 and A 0 i is the Jacobian of r 0 i with respect to w 0 i. The sub-iteration 
updates are applied with a solution dependent relaxation factor rjj = m/(Awj m + 2m). 

<oV = <oi + (13) 

At the beginning of the sub-iterations to = 1 we initialize woi with the average value of In pj at the cell 
centers bounding the surface. This initialization and use of relaxation factors r/j are designed to accommodate 
difficult transients as the global solution evolves. If sub-iterations fail to converge after 300 steps a warning 
is printed but the global solution is allowed to continue. It is not uncommon to see such warnings early in 
the evolution of the simulation but full convergence at all surface nodes is achieved as the global solution 
converges. 

2. Ablation with T w < T sw 

There is one extra unknown, to, when ablation is engaged. In contrast to the previous implementation the 
wall temperature in this case is computed after the coupled solution of Eqs. 1, 2, 3, and 9. 

A™ Aw 02 = r™ (14) 

rg?> is the same as roi but with the addition of one extra element equal to the residual of Eq. 9. In like 
manner, W 02 is the same as wqi but with the addition of one extra unknown m. A^ is the Jacobian of ro 2 
with respect to w 0 2 - The sub-iterations are engaged exactly as defined in the previous section so that 

O!) 1 = <02 + Vj Awj,o 2 1 <j<N s (15) 

< bt 1 = <02 + AtOj -,02 j = N s + 1 (16) 

There is no additional relaxation factor required on ?h for the sub-iterations to converge. However, at the 
conclusion of the sub-iterations we take a fraction of the updated blowing rate. 

to" + 1 = fAry (A r a+1 ) i02 + to" (17) 

The temperature is next updated using Eq. 8. 


rj-in -\- 1 


Qc+^ra d -^ n+1 K\ 1/4 


ea 


+ (1 - e)TZ 


(18) 


3. Ablation with T w > T sw 

In this case T w is computed from the requirement of chemical equilibrium with the surface using Eq. 9. 


K ch ar, C (T™) 


(Pc.oPcp 

M c 


1/2 


(19) 
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A sub-iteration over index to on this equation is executed using a Newton method until convergence is 
achieved. A fraction of this updated value is applied to advance the surface temperature. 


rjin+1 rrin , rpm 

^ w ^ w ' t'y-L w 


T n ) 

w ) 


(20) 


Densities are computed by the simultaneous solution of Eqs. 1, 2, and 3 for lnp^o with T w = T™ +1 and 
to = m" using the algorithm described previously in Eqs. 12 and 13. Finally, to" + 1 is computed from the 
surface energy balance Eq. 8. 


TO 


n+1 


= £ 


(<£ + - ^CC +1 ) 4 ) 


L»+ 1 


+ (1 — e)m n 


(21) 


B. Implementation of Surface Recession: Shape Change and Grid Adaptation 

1. Charring Ablator 

The charring ablator model defines mass loss in two ways. (Mechanical erosion is ignored.) As heat is con- 
ducted into virgin material a component of the material, usually a resin, undergoes decomposition through 
pyrolysis. If the back side of the TPS is sealed, pyrolysed gases escape into the shock layer through the 
remaining char. It is assumed that the original contour of the TPS is preserved until the char itself, predom- 
inantly carbon, begins to sublimate - a process that requires higher temperatures then pyrolysis. Changes to 
the outer mold line between two trajectory points are therefore a function of the char mass loss rate, rn c har , 
the char density p c har and the time between trajectory points, At. 

An = ™ charAt = ++_ At = ™ A t (22) 

Pchar Pv Pchar Pv 

The last equalities are enabled by the steady-state ablation approximation of Eq. 5. 


2. Non-charring Ablator 

In this model, only a single component ablation process is considered without any identifiable char remaining 
as surface recession proceeds. In the model setup, it is convenient to define a “char” density, p c har equal 
to the virgin TPS density, p v . Consequently, “char” refers to the receding material which, in this case, is 
identical to the virgin material. Therefore, from Eq. 5 we note that m g = 0 and therefore 


TO c ^ ar At 

Pchar 


m m . 

At = — At 


Pchar 


(23) 


This result is equivalent to the expression for the charring ablator with the additional simplification that 

Pchar / Pv — E 


3. Explicit Shape Change 

Blowing rates rh l: j are evaluated at cell wall i,j in the cell-centered, finite- volume LAURA code. The 
displacement at node (i+1/2, j + 1/2) is evaluated as the average of surrounding surfaces using the converged 
solution for trajectory point nt. 


/\n nt 

L * n i+l/2, j + 1/2 


t nt + 1 - t r 
4 p v 


« 


m, 


^nt 


m 


^nt 

i,j + 1 


m, 


^nt \ 

i+l,j+lj 


(24) 


The displacement components at this node are computed from the cross-product of the diagonals defined by 
the surrounding cells. 


r!r nt 

l,i+l/2, j+1/2 
dr nt 

“'2, i+1/2, j+1/2 
'h+1/2, j+1/2 

\r nt 

^ i+1/2, j+1/2 


r i+3/2,j-l/2 r i-l/2,j+3/2 


'i+3/2, j + 1/2 


_ ffn, t 


1/2, j — 1/2 

zfllt 


i+1/2, j+1/2 X ^^2, i+1/2, j+l/ 2 /W^l, i+1/2, j+1/2 X ^2, i+1/2, j+1/2 I 

^ n i+l/2,j+l/2^i+l/2,j+l/2 


(25) 

(26) 

(27) 

(28) 
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The volume grid for the new solution at time t nt+1 is now computed as 


rfllt+l _ zfilt , \znt 

, i+l/2,j'+l/2,fc+l/2 — 'i+ l/2,i + l/2,fe+l/2 "T ^'»+l/2,j+l/2 


(29) 


where each volume node is displaced according to its corresponding surface node in this structured grid 
formulation. In the case of axi-symmetric geometries displacements are computed along the original bounding 
planes defining the symmetry boundaries. 

Explicit shape change formulations are known to be subject to numerical instabilities that can induce 
oscillations in surface contours and heating . 2 These problems are most robustly dealt with through the use 
of a sub-iterative procedure that computes the shape at time nt + 1 as a function of the simulated solution at 
time nt + 1 . 2 Root causes of the explicit stability limit are explored in the Appendix and a filtering algorithm 
is introduced that enhances explicit stability. The baseline explicit update described above is referenced as 
a two-cell stencil algorithm. The enhanced explicit update described in the Appendix is referenced as a 
four-cell stencil algorithm. 


4- Volume Grid Adaptation 

All simulations incorporate grid adaptation as an integral part of the flow solver. The outer boundary of 
the grid is aligned with the captured shock. The near wall resolution requires 0.1 < pc ^ n < 1. In cases with 
shape change, the volume grid is re-adapted every 500 relaxation steps for the first 5000 iterations giving a 
total of 10 adaptations per trajectory point. 


V. Simulations 

All of the flow field simulations presented here, including fully coupled radiation and ablation, are executed 
with programs LAURA 16 and HARA. 1 ' The test cases span conditions for which peak mass loss rate ranges 
from 0.06 kg/m 2 -s to 140 kg/m 2 -s. Line implicit relaxation is applied for all solutions. 

A. IRV-2 

The IRV-2 geometry is described in the references . 2 The IRV-2 vehicle is a sphere-biconic-cylinder with a 
nose radius of 0.01905 m and total length of 1.386 m. The biconic angles are 8.42 deg. and 6.10 deg. The 
fore-cone breaks at 0.1488 m which marks the extent of the present simulations. Free stream conditions 
across the trajectory are defined in Table 1. A trajectory data file enables LAURA to sequentially compute 
a flow field at a trajectory point, compute a shape change based on the m distribution, re-grid, and move 
on to the next trajectory point. The baseline grid uses 64 cells normal to the body and 64 cells around the 
body with 30 of those cells on the nosetip. The outer boundary is aligned with the captured bow shock. 
The cell Reynolds number at the wall is set to 1.0. A grid convergence test (double the grid density across 
the shock layer) at trajectory point 16 indicates stagnation point heating converged to 0.75%. Nineteen 
species for products of high temperature air and products from a carbon ablator include: N, O, N%, O 2 , 
NO, N + , 0 + , N£ , O 2 , NO + , e~ , C, C 2 , C 3 , C 4 , C 5 , CO, CN, Laminar flow in thermochemical 

nonequilibrium is specified for the first 15 trajectory points. Thermal equilibrium is specified for trajectory 
point 16 and beyond. These specifications are consistent with the original refernce . 2 Two thermal protection 
system (TPS) densities are simulated for this trajectory: a low-density carbon ( p v = 544.63 kg/m 3 ) and 
carbon-carbon (p v = 1800 kg/m 3 ). The low-density tests induce larger recessions for equivalent blowing 
rates and provide greater opportunity to explore stability of the shape change algorithm. Carbon-carbon is 
the material specified for the IRV-2 test. 

1. Low-Density TPS 

Results of simulations for the first six trajectory points are presented in Fig. 1 for the low-density TPS. Mass 
loss rate rh is shown in cyan, the atomic carbon mass fraction cc is shown in green, the surface temperature 
is shown in blue, convective heating is shown in red, and body shape is shown in black as a function of time 
during first 11.5 seconds of IRV-2 test. The axes are defined to provide an aspect ratio of approximately one 
for the body shape (black line). The body shape is presented nose down. Recession of the stagnation point 
with time can be read as the point where the black line intersects the left axis. The carbon mass fraction is 
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tracked because its equilibration with the solid carbon on the surface is a key element of closing the surface 
boundary condition. 


Table 1. Free Stream Conditions for IRV-2 2 


Traj- 

Point 

Time 

00 

Altitude 

(m) 

Velocity 

(m/s) 

Temp. 

(K) 

Pressure 

(Pa) 

Density 

(kg/m 3 ) 

Mach 

No. 

Thermal 

State 

1 

0 

66935 

6780.6 

227.81 

8.1757 

1.2505xl0" 4 

22.41 

TNEQ a 

2 

4.25 

55842 

6788.3 

258.02 

37.362 

5.0454xl0 -4 

21.08 

TNEQ 

3 

6.75 

49290 

6785.2 

270.65 

88.118 

1.1344xl0" 3 

20.57 

TNEQ 

4 

8.75 

44042 

6773.0 

261.40 

169.50 

2.2593xl0" 3 

20.90 

TNEQ 

5 

10.25 

40108 

6752.4 

250.35 

287.14 

3.9957xl0 -3 

21.29 

TNEQ 

6 

11.50 

36836 

6722.0 

241.50 

445.52 

6.4268xl0 -3 

21.58 

TNEQ 

7 

12.50 

34229 

6684.3 

234.30 

644.52 

9.5832xl0" 3 

21.78 

TNEQ 

8 

13.25 

32283 

6644.9 

228.76 

863.14 

1.3145xl0" 2 

21.91 

TNEQ 

9 

13.95 

30480 

6596.7 

226.91 

1127.6 

1.7313xl0 -2 

21.84 

TNEQ 

10 

14.75 

28236 

6527.1 

224.73 

1568.1 

2.4310xl0 -2 

21.71 

TNEQ 

11 

15.50 

25772 

6428.3 

222.35 

2256.0 

3.5348xl0 -2 

21.50 

TNEQ 

12 

16.25 

22949 

6286.6 

219.47 

3520.8 

5.5888xl0 -2 

21.17 

TNEQ 

13 

17.00 

19790 

6091.7 

216.65 

5705.3 

9.1741xl0" 2 

20.64 

TNEQ 

14 

17.75 

16355 

5836.4 

216.65 

9723.1 

1.5635xl0 _1 

19.77 

TNEQ 

15 

18.25 

13962 

5631.8 

216.65 

14170 

2.2786xl0 _1 

19.08 

TNEQ 

16 

18.50 

12748 

5519.6 

216.65 

17379 

2.7946xl0 _1 

18.70 

TEQ fc 

17 

18.75 

11528 

5401.2 

216.65 

20984 

3.3743xl0 _1 

18.30 

TEQ 

18 

19.00 

10309 

5277.1 

221.31 

25309 

3.9840xl0 _1 

17.69 

TEQ 

19 

19.50 

7892 

5014.3 

236.86 

36168 

5.3196xl0 _1 

16.25 

TEQ 

20 

20.00 

5536 

4736.5 

252.11 

50198 

6.9366xl0 _1 

14.88 

TEQ 

21 

20.50 

3273 

4449.6 

267.04 

67922 

8.8610xl0 _1 

13.58 

TEQ 

22 

21.00 

1129 

4159.7 

280.68 

88253 

1.0954x10° 

12.38 

TEQ 

23 

21.28 

0 

4000.0 

288.15 

101325 

1.2250x10° 

11.75 

TEQ 


“ thermal nonequilibrium, b thermal equilibrium 

Fig. 2 compares simulation results at trajectory points 7, 9, and 11 using the baseline, 2-cell stencil, 
explicit shape change algorithm (Eqs. 24 - 27) with the 4-cell stencil algorithm (Eq. 32 defined in the 
Appendix). The two-cell stencil results are on the left and the four-cell stencil on the right for a given 
trajectory point. Oscillations in heating and blowing rate are significantly larger for the baseline algorithm. 
Oscillations in body shape are present but are difficult to see except at the 15.5 s point for the baseline 
algorithm. The correlation of heating spikes with local surface compressions is more evident in Fig. 3. This 
figure shows the evolution of the node shape (black) and heating (red) over the first eleven trajectory points 
with the low-density (TPS). Very slight surface waves produce a relatively large perturbation in surface 
heating, blowing rate, and atomic Carbon mass fraction. Surface temperature is relatively insensitive to 
these perturbations. Averaging over a larger stencil inhibits the growth of instability but does not eliminate 
it. In essence, a larger effective Ax yields a larger allowable At as discussed in the Appendix. 

2. Carbon-carbon TPS 

Results of simulations for odd number trajectory points using the carbon-carbon TPS are presented in Fig. 
4. Note that the first three subfigures use a different range than the last nine subfigures. Shape change for 
these cases used the 4-cell stencil algorithm. Note that the mass loss rate m at a given trajectory point for 
the high-density TPS (Fig. 4 d, e, f) is approximately equal to the mass loss rate for the low-density TPS 
(Fig. 2). Even though the evolution of m is equivalent up to trajectory point 11 at 15.5 s for the two TPS 
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(a) t = 0 s (b) t = 4.25 s 




(c) t = 6.75 s 


(d) t = 8.75 s 




(e) t = 10.25 s (f) t = 11.50 s 

Figure 1. Mass loss rate m (cyan), atomic carbon mass fraction cc (green), surface temperature (blue), 
convective heating (red), and body shape (black) as function of time during first 11.5 seconds of IRV-2 low 
density char test. 
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(a) t = 12.5 s, 2-cell stencil 



(b) t = 12.5 s, 4-cell stencil 




(c) t = 13.95 s, 2-cell stencil (d) t = 13.95 s, 4-cell stencil 




(e) t = 15.5 s, 2-cell stencil (f) t = 15.5 s, 4-cell stencil 

Figure 2. Mass loss rate rh (cyan), atomic carbon mass fraction cc (green), surface temperature (blue), 
convective heating (red), and body shape (black) as function of time for 12.5 < t < 15.5 IRV-2, low-density TPS. 
Slight irregularities in the body shape cause large oscillations in heating that are moderated through use of a 
filtering algorithm extending over a larger stencil. 
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Figure 3. Evolution of body shape and heating distribution for first 11 simulation points spanning 0 < t < 15.5s 
with low-density TPS. 


systems one observes that the high-density TPS shows no sign of a shape change induced instability (Fig. 
4(f)) whereas the the low-density TPS case shows onset of the instabilities at this same point (Fig. 2(f)). 
This behavior indicates that the shape change instability is associated with the recession rate distribution as 
opposed to mass loss rate distribution. Some weak oscillations in the heating and blowing rate are evident 
at trajectory point 17 (18.75 s, Fig. 4(i)). These oscillations grow slowly through trajectory point 19 (t = 
19.5 s, Fig. 4(j)). Significant oscillations first appear at trajectory point 20 (not shown) and continue to 
trajectory point 23 (Fig. 4(1)). 

A key result of Fig. 4 is that it demonstrates that boundary condition implementation is robust for a 
wide range of conditions. Temperature and heating rate vary from 1472 K and 24 W/cm 2 at the trailing 
edge for trajectory point 1 to 4450 K and 5000 W/cm 2 at the stagnation point for trajectory point 17. 
Trajectory points 1-13 have the switch temperature T sw = 2800K included in the domain. In these cases, 
different relaxation sequences are used to compute the boundary conditions. 

The convergence history for this sequence of cases is presented in Fig. 5. Trajectory point 1 had been 
computed previously. Trajectory points 2-15 (Fig. 5(a)) were run assuming thermal non-equilibrium. 
Trajectory points 16 - 23 (Fig. 5(b)) were run assuming thermal equilibrium. The cases are automatically 
sequenced using a trajectory data file. A detail of the convergence history for trajectory points 8 and 9 is 
presented in Fig. 5 c. A spike in the error norm is evident after sequencing to the next trajectory point. Ten 
additional spikes associated with grid adaptation every 500 relaxation steps are also recorded. There is a 
nearly monotone decrease in error norm after the last grid adaptation until a convergence criteria of 3 10“ 10 
is achieved. (The final two trajectory points show the error norm reduction stalls - probably associated with 
the growing shape change instabilities.) The non-monotone spikes at approximately 39500 s and 44000 s 
are associated with a change of temperature across T sw late in the convergence process causing an abrupt 
change in the relaxation of the boundary condition. (The boundary conditions do not change - only the 
relaxation process to satisfy the boundary conditions change.) Still, one can see in the detailed Fig. 6 that 
all surface quantities are smoothly resolved across T sw = 2800K at convergence. Note that cc rapidly goes 
to trace levels at T = 2800K and that to remains at a small, nearly constant level for T w < 2800K. 

At trajectory points 6-9 a negative value of to is computed in regions just behind the blunt nose. Surface 
pressure and temperature tend to rapidly decrease in these regions. A significant amount of ablated species 
are carried to these points from upstream locations. Carbon may condense on the surface (negative to) in 
order to equilibrate the atomic carbon vapor pressure. 18 

Stagnation point results are compared to a reference calculation by Kuntz, et. al. 2 in Fig. 7. The 
LAURA results in red with diamond symbols are overlayed on the original figures pulled from the reference. 
The relevant reference curves are indicated by solid black lines with filled square symbols indicating their 
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-1 4000 0.008 



(a) t = 0.00 s, point 1 



(b) t = 6.75 s, point 3 



(c) t = 10.25 s, point 5 





(d) t = 12.50 s, point 7 


(e) t = 13.95 s, point 9 


(f) t = 15.50 s, point 11 




(g) t = 17.00 s, point 13 


(h) t = 18.25 s, point 15 


(i) t = 18.75 s, point 17 




(k) t = 20.50 s, point 21 



(1) t = 21.28 s, point 23 


Figure 4. Mass loss rate m 
convective heating (red), and 


(cyan), atomic Carbon mass fraction cc (green), surface temperature 
body shape (black) as function of time for carbon-carbon TPS. 


(blue), 


13 of 23 


American Institute of Aeronautics and Astronautics 


kg/m 3 - 




(a) trajectory points 2 - 15, TNEQ (b) trajectory points 16 - 23, TEQ 



(c) trajectory points 8 and 9 


Figure 5. Convergence history. 




(a) t = 13.25 s, point 8 (b) t = 13.95 s, point 9 

Figure 6. Detail of boundary condition distributions highlighting behavior across T w = T sw = 2800K. 
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iterative algorithm. Solid black circles indicate results from the engineering code ASCC 19 used by Kuntz 
et. al. as another reference calculation. The reference calculation, SACCARA, used a material response 
code, COYOTE, with initial conditions corresponding to a cold, non-ablating wall. The LAURA simulation 
replaces the material response code with the steady state ablation approximation - a simpler but less accurate 
model of the true evolution of heat transfer through the TPS and subsequent mass loss. Note that an earlier 
LAURA simulation 20 of this case did show good agreement with the early material response. That simulation 
coupled a one-dimensional material response code and was not used here because it implicitly incorporates 
approximations that bypass the boundary condition relaxation strategy tested here. 

The LAURA surface temperature starts at a very high temperature because the body is initialized at 
time t = 0 as if it had already achieved a steady state ablation condition. (Fig. 7(b) ) It takes the more 
accurate material response simulation approximately 10 seconds before the surface temperature asymptotes 
to the level predicted by the steady state ablation approximation. The surface temperature is computed 
from Eq. 9 requiring equilibration of atomic carbon with surface char. 

The LAURA heating rate starts 17% lower than the SACCARA value at t = 0 but after 10 seconds 
begins to rise above the reference calculation. (Fig. 7(a) ) This difference is again associated with the 
LAURA initial condition that has a fully established blowing rate and larger surface temperature at t = 0 as 
compared to the reference. Peak heating occurs near t = 19s for both LAURA and SACCARA with LAURA 
approximately 20% larger than the reference value. The blowing rate predicted by LAURA is larger than 
that predicted by SACCARA/COYOTE up to the onset of peak heating at approximately t = 18s. (Fig. 
7(c) ) It is likely that the steady state ablation approximation is still too aggressive during this period and 
the conduction energy flux into the TPS is not yet balanced by the energy flux in the outgassing. The higher 
heating rates in LAURA compared to SACCARA beyond t = 12s are more difficult to explain. While there 
has been more surface recession in LAURA at a given time the body shape appears to be slightly blunter 
than SACCARA predicts so effective radius is not a satisfactory explanation. Both simulations compute 
an equilibrium species mass fraction at the surface. A grid convergence test at trajectory point 16 using 
128 cells across the shock layer produced a 0.76% reduction in heating - not nearly enough to account for 
current differences. LAURA was run with ionization but the mass fraction of NO + in the stagnation region 
is only 0.0004. At present it is thought that the conduction terms and energy storage terms ignored in the 
steady state ablation approximation remain sufficiently large compared to energy exiting the surface through 
sublimation for this simulation. 
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(a) heating 


(b) surface temperature 




(c) mass loss rate (d) recession 

Figure 7. Stagnation point heat flux, temperature, blowing rate, and recession. Black lines taken from Kuntz, 
et. al. for material response using COYOTE coupled with the SACCARA flow solver. Red lines are for 
LAURA using quasi-steady, equilibrium ablation model. 
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B. Galileo Probe 


The Galileo probe entered the Jovian atmosphere at a relative velocity of 47.4 km/s in December, 1995. 21 
The vehicle forebody shape is a 44.86 deg. spherically capped cone with nose radius 22.2 cm and base radius 
of 63.2 cm. Free stream conditions for the current simulation are: V/o = 41591 m/s, p^ = 3.49ccl0 -4 kg/m 3 , 
and Tqo = 300. K. The flow field is resolved with 52 cells around the body and 100 cells across the shock 
layer. Sixteen species for products of high temperature hydrogen-helium atmosphere and ablation products 
include: H 2 , H, H + , He , He + , e~ , C, C + , CO, C 2 , C 3 , C 2 H, O, 0 + , C 5 , C0 2 . The elemental mass 
fractions of a charring ablator are specified: cc = 0.92, ch = 0.02, and co = 0.06. The free stream 
atmospheric composition is specified c # 2 = 0.76056 and Cn e = 0.23944. The shock layer is assumed to be 
in chemical equilibrium using the free-energy-minimization option in LAURA. The Cebeci-Smith algebraic 
turbulence model is applied. Radiation is fully coupled with the solution of the shock layer. 

Figure 8 presents distributions of mass loss rate rh in cyan, the atomic Carbon mass fraction cc in 
green, the surface temperature in blue, convective heating in red, and radiative heating in violet. The 
radiative heating is much larger than convective in this case. Peak mass loss rate exceeds 7 kg/m 2 -s. The 
surface temperature is nearly flat at approximately 4000 K. The solid line indicates use of a film coefficient 
approximation which is in excellent agreement with the new boundary condition implementation shown by 
dashed lines. 

The large blowing rates driven by radiative heating essentially displace the boundary layer off the body. 
The extent of this displacement is evident in Fig. 9(a) which shows the carbon mass fraction in the blown 
gas layer and in Fig. 9(b) which shows the ionized atomic hydrogen in the radiating, inviscid layer. 



Figure 8. Mass loss rate rh (cyan), atomic Carbon mass fraction cq (green), surface temperature (blue), con- 
vective heating (red), and radiative heating (purple) at 51.16 s (peak heating) for the Galileo probe simulation. 
The use of film coefficient approximation (solid line) is in excellent agreement with the new boundary condition 
implementation (dashed line). 


C. 10 km/s at Ground Level 

Conceptual ideas for launching small satellites from the ground using electromagnetic ring accelerators require 
a thermal protection system for the carrier vehicle projectile (CVP) capable of withstanding the aerothermal 
environment associated with a velocity of approximately 10 km/s at ground level densities . 22,23 Exploratory 
simulations are executed here to shake out any numerical issues associated with the large ablation rates 
computed for these environments. Similar studies have been engaged for related concepts . 24 The main 
result of this section is that the boundary relaxation algorithm worked without problem in this challenging 
simulation. 

The vehicle proposed here, derived from earlier concepts , 25,26 is a 7 deg. spherically capped cone with 
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(a) C 



(b) H+ 


Figure 9. Contours of atomic Carbon and ionized atomic Hydrogen indicating extent of the blown gas layer 
and inviscid layer including streamlines from the free stream and from the body for Galileo. 


5 cm nose radius and 3 meter length. Free stream conditions at ground level are: V^o — 10000 m/s, 
Poo = 1.225 kg/m 3 , and = 288.1 K. The flow field is resolved with 60 cells around the body and 64 cells 
across the shock layer. Twenty-six species for products of high temperature air and products from an Avcoat 
ablator include: N, O , N%, O 2 , NO, N + , 0 + , N.f, , NO + , e~, C, C + , CO, CO 2 , C 2 , C 3 , C 2 H, CN, H, 
H + , H 2 , C 2 H 2 , C 5 , HCN , CH. The elemental mass fractions of pyrolysis gases are specified: cc = 0.547, 
ch = 0.093, co = 0.341, and cn = 0.019. The elemental mass fraction of the char is cc = 1.0. The char 
density is 256.3 kg/m 3 . The virgin ablator density is 544.6 kg/m 3 . (We make no assertion that this is the 
right ablator for this task or that steady-state ablation is an appropriate approximation - these are simply 
starting points for this study.) Both thermochemical equilibrium and non-equilibrium models have been 
tested. Only the equilibrium results are presented with coupled radiation and turbulent boundary layer. 




Figure 10. Comparison of aerothermal loads with and without ablation for the carrier vehicle projectile (CVP). 


The severity of 10 km/s at ground level should not be understated. Consider the environments described 
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200000 


200 -.5000 



Figure 11. Mass loss rate m (black), surface temperature (blue), convective heating (red), and radiative heating 
(violet) for the carrier vehicle projectile (CVP) using two different turbulence models. 


in Fig. 10 which highlights the changes to heating and shear associated with ablation. First note that 
pressure distributions are insensitive to ablation. The stagnation pressure in Fig. 10(b) at 116 million N/m 2 
(116 MPa, 16824 lbf/in 2 ) is equivalent to the pressure at a depth of over 7 miles below sea level in the ocean. 
The pressure on the seven degree conical section rapidly drops to approximately 2 million N/m 2 (2 MPa, 
290 lbf/in 2 ) which is equivalent to the pressure at a depth of 0.1 miles below sea level in the ocean. Shear 
stresses drop by approximately a factor of two due to the presence of ablation. The peak shear occurs near 
the sphere cone junction. Peak radiative heating levels occur at the stagnation point (Fig. 10(a)). The 
Cebeci-Smith algebraic turbulence model is applied everywhere (transition begins at the stagnation point) 
and so the peak convective heating levels occur just ahead of the sphere cone junction. The peak heating 
levels are more than a factor of 10 higher than the peak radiative heating encountered in the previous Galileo 
simulation (Fig. 8). Note that convective heating drops by an approximate factor of two due to ablation. The 
radiative heating with ablation is less than the non-ablating case in the stagnation region but is larger than 
the non-ablating case over the cone flank. This increase is associated with the presence of high temperature 
ablation products in the shock layer. These higher radiative heating levels are still almost a factor 10 less 
than the corresponding convective heating levels on the flank. The flank total heating levels are slightly less 
than the peak Galileo levels in Fig. 8. However, the Galileo pressure level at the peak heating location (not 
shown) is approximately a factor of 4 lower than the CVP pressure level on the flank. 

Fig. 11 shows mass loss rate to in black, the surface temperature in blue, convective heating in red 
and radiative heating in violet for two turbulence models. There is no basis to assign uncertainty to any 
of these results; the two turbulence models are tested to provide an initial assessment of sensitivity. The 
baseline algebraic Cebeci-Smith model is shown with solid lines and Wilcox’s k-omega model (2006) with 
compressibility correction is shown with dashed lines. The stagnation point heating is 45000 to 65000 W/cm 2 
and the associated blowing rate is 75 to 95 kg/m 2 -s, depending on the turbulence model. The largest 
difference between the models is just ahead of the sphere-cone junction where the Cebeci-Smith model is at 
125000 W/cm 2 and the k-omega model is at 170000 W/cm 2 . The associated peak blowing rate varies from 
approximately 125 to 145 kg/m 2 -s. A shape change computation has not been attempted for this case but it 
is interesting to note that an 1800 kg/m 3 carbon-carbon ablator would experience an 8 cm/s recession rate 
- not counting any mechanical erosion. The k-omega results are thought to be too high because the model 
admits some production of turbulent kinetic energy behind the normal shock in this case. The stagnation 
point radiative heating levels are 103000 to 115000 W/cm 2 . The surface temperature is relatively insensitive 
to the turbulence model and varies from approximately 4800 K at the stagnation point to 4100 K on the 
flank. 

While the magnitude of the blowing rate is very large in this example the peak dimensionless value of 
m/ (pooVoo) is only 0.0118 which is approximately the stagnation point value of the IRV2 at trajectory point 
8. Consequently, the blown gas layer sits very close to the body. The blowing rate is high but the pressure 
is also high which limits the extent of the ablated gases. For example, the shock standoff distance over the 
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50 mm radius body is 15.12 mm at the sphere-cone junction. (Cebeci-Smith turbulence model) The dividing 
streamline separating free-stream gases from ablated gases is 0.42 mm above the surface. There are 34 cells 
stretched across this inner portion of the boundary layer. The mesh size at the wall is 1.16 10~ 4 mm. 

The design of a CVP will require a more comprehensive coupled material response - one that includes 
effects of temperature on compressive strength. Like the earlier IRV2 simulations, the vehicle is not “born” 
with a steady state ablation condition. The initial heating and shear loads will not be moderated by an 
ablating boundary condition. As altitude increases, free-stream density and velocity decrease and so the 
most severe conditions persist only on the order of seconds. A nose tip must be formulated that will be 
consumed by the flight environment in a predictable manner, lasting just long enough to protect the CVP 
flank and payload until it exits the sensible atmosphere. The flank of the CVP also sees severe aerothermal 
loads but these loads are not too dissimilar from those encountered in the Galileo probe. These engineering 
tasks are challenging but no show-stoppers are evident at this early stage. 

VI. Concluding Remarks 

A relaxation algorithm to implement equilibrium, steady-state ablation boundary conditions is defined 
and tested for the purpose of providing strong coupling with a hypersonic flow solver. The objective is to 
remove correction factors or film cooling approximations that are usually applied in uncoupled or loosely 
coupled implementations of the flow solver and the ablation response. Three test cases are considered - the 
IRV-2, the Galileo probe, and a notional slender, blunted cone (carrier vehicle projectile, CVP) launched at 
10 km/s from the Earth’s surface. The algorithm recognizes that critical elements of the computed bound- 
ary condition, the vapor pressure of carbon at the surface and the magnitude of the blowing rate, can vary 
by orders of magnitude as a function of surface temperature. Robust convergence of the linearized set of 
boundary conditions in this circumstance may be compromised. Consequently, a successive substitution is 
employed and the order of succession is varied as a function of surface temperature to obtain converged 
solutions. The algorithm is tested under the assumption of equilibrium, steady-state ablation. This approx- 
imation introduces some simplifications into the surface energy balance that obviates a need to engage a 
material response code. In future work, coupling with material response code requires minor changes to the 
surface energy balance equation and no significant change to the present boundary condition relaxation is 
anticipated. In like manner for future work, non-equilibrium ablation conditions will require replacement of 
partial equilibrium relations with species continuity in the boundary condition equation set. 

The implementation is tested on a specified trajectory for the IRV-2 to compute shape change under 
the approximation of steady-state ablation and one-dimensional material response with a shock layer in 
thermochemical non-equilibrium. This trajectory involved surface temperature variations that tested all 
options in the temperature dependent relaxation algorithm. A convergence history for each point of the 
trajectory is presented. Convergence for each point was achieved in 15000 to 30000 relaxation steps, including 
10 grid adaptations during the first 5000 steps using a boundary condition relaxation factor e = 0.001. Two 
TPS densities were tested to investigate issues associated with instabilities arising from shape change. An 
approximate stability analysis indicates that instabilities arise from the dependence of the recession rate on 
local perturbations to the body surface contour as opposed to absolute levels of recession rate. A four-cell 
stencil is defined to enhance stability of the explicit shape change algorithm. Comparisons are made to a 
previous study which shows, not surprisingly, that a steady-state ablation approximation over-predicts the 
recession rate with largest differences encountered early in the trajectory. 

A simulation on the Galileo probe, in which the heating rate is almost entirely due to radiation with 
the boundary layer completely blown off the body, was implemented with 16 species in thermochemical 
equilibrium. Fully coupled radiation and an algebraic Cebeci-Smith turbulence model were engaged for 
this case. Comparisons to an earlier result using film coefficient B-prime approximation showed excellent 
agreement with the new relaxation algorithm. 

The last simulation addressed a most extreme case, 10 km/s at sea level, to test the boundary relaxation 
algorithm at very high blowing rates. The case describes a notional concept for a carrier vehicle projectile 
(CVP) flung into orbit from an electromagnetic ring launcher. Here again, the solution exhibits good 
convergence with specifications of an Avcoat ablator including 26 species in thermochemical equilibrium, 
radiation, and a coupled algebraic or two-equation turbulence model. Stagnation pressures on the order of 
116 MPa and peak heating rates in excess of 125000 W/cm 2 are computed for this case. 


20 of 23 


American Institute of Aeronautics and Astronautics 



Appendix - Shape Change Stability Analysis 


Inspection of Fig. 3 suggests a correlation between the appearance of compression surfaces and localized 
hot spots. If one passes a reference plane through any subsection of the ablating nose using least squares 
fitting with perpendicular offsets 27 one notes that, relative to this reference plane, valleys and falling surfaces 
(expansions) tend to exhibit lower heating while peaks and rising surfaces (compressions) tend to exhibit 
higher heating. The correlation is imperfect and degrades as features grow and become more jagged - 
engaging non-linear interactions more strongly or feeling effects of upstream shock curvature. Still, the 
initial behavior suggests a simple analytic model for evaluating the stability of shape change. 

Consider two-dimensional flow over an ablating surface. A plane is fit through a segment of the surface 
and surface deformations are defined locally by y(x) where y is measured relative to the plane and x varies 
between 0 and 2tt as shown in Fig. 12. Assume that the recession rate y can be expressed as a function 
of an average recession rate over the segment, yo, with perturbations to this average average recession rate 
defined as a function of local slope and curvature. Thus, 


y = a 


dx 2 


(A 

ax 


y o 


(30) 



Figure 12. Response of notional, sinusoidal surface perturbation to ablation. Black lines and symbols refer 
to original surface and surface nodes at t = 0. Red lines and symbols refer to receded surface and nodes at 


t = 0.4935. 


The analytic solution to Eq.30 is given by 

N 

y{x , t) = yot + ^ e~ an * cos(n(x + bt)) + B n sin(n(x + bt))^j (31) 

n = 0 

and the initial shape y(x,0) may be defined with a Fourier series. The requisite assumption of periodic 
continuation of surface deformations to adjacent segments with trigonometric basis functions is unrealistic 
but the simplification enables some insights into stability of the explicit shape change algorithm. As an 
example, consider the solution to Eq. 30 with a = 2, b = 3, yo = ^3 with y(x,0) = sin(x), the black line 
in Fig. 12. The red line indicates the solution at t = 0.4935. In this example, the features of the black 
curve have moved to the left a distance bt, the wave amplitude has decreased by a factor e~ at , and average 
level of the curve has receded to y = yot. If a < 0 on a segment then the the recession rate in valleys is 
larger than the recession rate on peaks; consequently, peaks are not ablated as quickly as valleys and the 
amplitude of perturbations grow. If a > 0 then the amplitude of surface perturbations goes to zero and the 
recession proceeds at a constant rate yo- A numerical solution of Eq.30 using backward-time, centered-space 
discretization is subject to the same stability limits on time step as the wave equation and heat conduction 
equation - At < min(4p, ^-). In this idealized problem the stability of an explicit shape change algorithm 
is therefore not a function of the average recession rate yo but instead is a function of the magnitude of 
perturbations of the recession rate relative to perturbations in the body shape. 
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A semi-analytic solution to Eq.30 can be generated using Eq. 31 that has no time step constraints. It 
exploits the fact that the parameters a, b , and yo are constant in time and it generates a Fourier series to 
recompute constants A n and B n for the perturbed shape after each time step. Preliminary attempts to 
extend this idea to evaluate shape change in LAURA have failed to produce a robust algorithm. The failures 
stem from the underlying assumption that the variation of localized heating rates can be described strictly 
as a function of local surface perturbations. One finds that least squares fits to recession rates on a segment 
can yield negative values for a and that the constants a, b , and j/o change rapidly and erratically from one 
segment to the next. While there may yet be a filtering algorithm to overcome these challenges some simple 
updates based on these ideas have served to enhance stability of the explicit shape change algorithm. 

The axisymmetric version of the shape change algorithm is defined here. Multi-dimensional extensions 
are theoretically obvious but have not yet been programmed or tested. Every node on the surface is updated 
as a function of two neighboring nodes to the right and left (five nodes along the surface defining four 
cell walls along the surface). Determine the equation of a reference line passing through these five nodes 
according to the least squares - perpendicular offsets algorithm. 2 ' Compute the area average blowing rate 
over the four cell faces. 


™' i + 1/2 — 

The recession at node i + 1/2 is now given by 

An i+1/2 = A tm i+1/2 /p v (33) 

The recession Anj_|_i /2 is made in a direction normal to the reference line described above. Evaluating 
displacements relative to these reference lines inhibits the formation of kinks in the surface and is consistent 
with the quasi-one-dimensional, steady-state ablation models used here. 



riyAj 



(32) 
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